IRIX Base Documentation 2002 November

home *** CD-ROM | disk | FTP | other *** search

/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / libblas / rotmg.z / rotmg

Wrap

Text File | 2002-10-03 | 15.8 KB | 342 lines

ROTMG(3F) Last changed: 11-2-98 NNAAMMEE SSRROOTTMMGG, DDRROOTTMMGG - Constructs a modified Givens plane rotation SSYYNNOOPPSSIISS Real CCAALLLL SSRROOTTMMGG ((_d ,, _d ,, _b ,, _b ,, _r_p_a_r_a_m)) _1 _2 _1 _2 Double precision CCAALLLL DDRROOTTMMGG ((_d ,, _d ,, _b ,, _b ,, _r_p_a_r_a_m)) _1 _2 _1 _2 IIMMPPLLEEMMEENNTTAATTIIOONN IRIX systems DDEESSCCRRIIPPTTIIOONN These routines compute the elements of a modified Givens plane rotation matrix. These routines have the following arguments: _d First diagonal element. (input and output) _1 SSRROOTTMMGG: Real. DDRROOTTMMGG: Double precision. On input, this value is the first diagonal element of the scaling matrix _D. On the first call to the routine, this value is typically 1.0. Subsequent calls typically use the value from the previous call. On output, this value is the first diagonal element of the updated scaling matrix _D''. _d First diagonal element. (input and output) _2 SSRROOTTMMGG: Real. DDRROOTTMMGG: Double precision. On input, this is the first diagonal element of the scaling matrix _D. On the first call to the routine, this value is typically 1.0. Subsequent calls typically use the value from the previous call. On output, this value is the first diagonal element of the updated scaling matrix _D''. _b x-coordinate. (input and output) _1 SSRROOTTMMGG: Real. DDRROOTTMMGG: Double precision. On input, this value is the _x-coordinate of the vector used to define the angle of rotation, before scaling (multiplying by the matrix _D). On output, this value is the _x-coordinate of the rotated vector, before scaling (multiplying by the matrix _D''). _b y-coordinate. (input) _2 SSRROOTTMMGG: Real. DDRROOTTMMGG: Double precision. On input, this value is the _y-coordinate of the vector used to define the angle of rotation, before scaling (multiplying by the matrix _D). It is unchanged on output. _r_p_a_r_a_m Real array of dimension 5. (output) SSRROOTTMMGG: Real array. DDRROOTTMMGG: Double precision array. This array contains rotation matrix information. The routine sets up the computed elements in _r_p_a_r_a_m from inputs _d , _d , _b , and _b . _1 _2 _1 _2 SSttaannddaarrdd GGiivveennss RRoottaattiioonn A standard Givens rotation (see SSRROOTTGG(3F)) is based on an orthogonal matrix _G that rotates points on a Cartesian _x_y-coordinate plane. To calculate the rotation matrix, you must provide the angle of rotation desired, or, equivalently, a vector (point) that lies along the angle of rotation. For a given planar point (_x , _y ), _G is formed so that: _r _r _ _ _ _ _ _ _ _ | x' | | c s | | x | G | x | | | | | | r | | r | | 0 | = |-s c | | y | = | y | | | | | | r | | r | - - - - - - - - 2 2 where _x' = sqrt (_x + _y ) . _r _r With this rotation matrix _G, you can then convert any number of existing planar points to the new (rotated) _x_y-coordinate system. For _n points, the rotations would be as follows: _ _ _ _ _ _ | x | | c s | | x | | i | | | | i | | 0 |<- |-s c | | y | | i | | | | i | - - - - - - for _i = 1, 2, ..., _n MMooddiiffiieedd GGiivveennss RRoottaattiioonn The algorithm for these routines is based on the following observation. The rotation matrix _G can be factored into a _s_c_a_l_i_n_g matrix (diagonal matrix) and _m_o_d_i_f_i_e_d rotation matrix _H, for which either the diagonal or the off-diagonal elements are units (that is, +1 or -1). Thus, to perform _m modified (scaled) rotations on _n planar points, requires only 2_n_m, rather than 4_n_m multiplications for the standard rotation. Because you may want to perform several successive rotations, this routine assumes that you have leftover scaling factors from your previous modified Givens rotation; that is, the routine requires you to input not only a planar rotation vector (_b , _b ) but also the squares of the diagonal elements of the scali_n1g m_a2trix, _d and _d . The actual rotation vector is specified as follows: _1 _2 _ _ _ _ _ _ 1/2_ _ | _x_r | | sqrt(_d_1) 0 | | _b_1 | _D | _b_1 | | _y_r | = | 0 sqrt(_d_2) | | _b_2 | = | _b_2 | - - - - - - - - where _d_1 and _d_2 are the input scaling factors. Given these inputs, the routine generates a new modified rotation matrix _H with units for either the diagonal or off-diagonal elements, and new elements _d_1' and _d_2' for the new scaling factors, and rotated but unscaled vector (_b_1', 0), with the following results: _ _ 1/2_ _ 1/2 _ _ 1/2_ _ _ _ _G | _x_r | _G _D | _b_1 | _D' _H | _b_1 | = _D' | _h_1_1 _h_1_2 | | _b_1 | | _y_r | = | _b_2 | = | _b_2 | | _h_2_1 _h_2_2 | | _b_2 | - - - - - - - - - - 1/2_ _ _ _ _D' | _b_1' | | _x' | = | 0 | = | 0 | - - - - where: 1/2 _D' ( = diag {sqrt(_d_1'), sqrt(_d_2')} ) uses the updated scaling factors _d '' and _d '', which are _d and _d on output. _1 _2 _1 _2 _H is stored in the output array argument _r_p_a_r_a_m; _b_1' is stored as _b_1 on output. 1/2 1/2 _D'' _H equals _G _D , nnoott _G, as implied earlier. You must account for the old scaling factors when calculating the new scaling factors. After calculating the matrix _H by using SSRROOTTMMGG/DDRROOTTMMGG, you can then use it in SSRROOTTMM/DDRROOTTMM to convert points to the new coordinate system. NNOOTTEESS MMeeaanniinngg ooff tthhee OOuuttppuutt VVaalluueess The output values are returned through arguments _d , _d , _b , and _r_p_a_r_a_m. _1 _2 _1 The scaling factors _d and _d are updated with each call to the routine. Although SSRR_OO1TTMM/DDRROO_TT2MM does not need the updated factors, they are needed in two other important contexts: * As input for subsequent calls to this routine. * As scaling factors for rotated but unscaled points (_x , _y ), which are output from SSRROOTTMM./DDRROOTTMM. _i _i In this second usage, the actual (scaled) points would be given by (sqrt(_d_1)*_x(_i), sqrt(_d_2)*_y(_i)). Doing this operation frequently on all your points is counterproductive. The main advantage of the modified rotation algorithm is to reduce the number of operations. If you fold in the scaling factors after each rotation, you are performing the same number of operations as in the standard Givens rotation. These two uses for the scaling factors are mutually exclusive; that is, if you fold the scaling factors back into all your points, you no longer need those factors for these routines. After folding in the scaling factors, and before the next call to SSRROOTTMMGG/DDRROOTTMMGG, reset _d and _d to 1.0. _1 _2 On output, _b represents the new _x-coordinate after rotating (but before scali_n1g) the rotation vector. Although the _y-coordinate of this vector is 0.0 (see the the previous discussion), the corresponding value _b is unchanged on output. _2 The output array argument _r_p_a_r_a_m specifies the format of matrix _H, and it holds the nonunit values of _H. This is the only output of these routines that SSRROOTTMM/DDRROOTTMM requires. Each element of _r_p_a_r_a_m has a specific meaning, as follows: _r_p_a_r_a_m(1) a flag parameter that specifies how the matrix is stored = 0.0 Off-diagonal elements of _H are units. = 1.0 Diagonal elements of _H are units. = -1.0 Rescaling case (see the following subsection). = -2.0 _H is the identity matrix; no rotation needed. _r_p_a_r_a_m(2) = _h if needed _1,_1 _r_p_a_r_a_m(3) = _h if needed _2,_1 _r_p_a_r_a_m(4) = _h if needed _1,_2 _r_p_a_r_a_m(5) = _h if needed _2,_2 CCaallccuullaattiinngg tthhee OOuuttppuutt VVaalluueess The following presents the algorithm for calculating the output values _d , _d , _b , and _r_p_a_r_a_m, based on the input values _d , _d , _b , and _b_2. T_h1is _a2lgo_r1ithm is presented without explanation or _p1roo_f2. _F1or a more complete discussion of the modified Givens algorithm, see the papers listed in the SEE ALSO section. CCaassee 11:: _b = 0 (trivial case) _2 In this case, no rotation is needed. The flag value, _r_p_a_r_a_m(1), is set to -2.0. When passed to routine SSRROOTTMM/DDRROOTTMM, this flag indicates that it should not do any rotations. On output from SSRROOTTMMGG/DDRROOTTMMGG, _d_1, _d_2, and _b_1 are unchanged. CCaassee 22:: | sqrt(_d_1)*_b_1 | > | sqrt(_d_2)*_b_2 | (|_x_r| > |_y_r|) In this case, the diagonal elements of _H (_h and _h ) are set to 1.0 . Thus, the _r_p_a_r_a_m values set on outpu_t1,_a1re as _f2,o_l2lows: _r_p_a_r_a_m(1) <-- 0.0 _r_p_a_r_a_m(3) <-- _h_2_1 = -_b_2/_b_1 _r_p_a_r_a_m(4) <-- _h_1,_2 = (_d_2*_b_2)/(_d_1*_b_1) The output values of the scaling factors are as follows: _d_1 <-- _d_1' = _d_1/_u _d_2 <-- _d_2' = _d_2/_u where _u = det(_H) = 1 + (_d_2*_b_2*_b_2)/(_d_1*_b_1*_b_1) . The output value of _b is as follows: _1 _b_1 <-- _b_1' = _b_1*_u If rescaling is needed, SSRROOTTMMGG/DDRROOTTMMGG will further modify these output values before the end of the routine. See case 4 later in this subsection. CCaassee 33:: | sqrt(_d_1)*_b_1 | <= | sqrt(_d_2)*_b_2 | (|_x_r| <= |_y_r|) In this case, the off-diagonal elements of _H are units (to be specific, _h_2,_1 = -1 and _h_1,_2 = 1). Thus, the _r_p_a_r_a_m values set on output are as follows: _r_p_a_r_a_m(1) <-- 1.0 _r_p_a_r_a_m(2) <-- _h_1,_1 = (_d_1*_b_1)/(_d_2*_b_2) _r_p_a_r_a_m(5) <-- _h_2,_2 = _b_1/_b_2 The output values of the scaling factors are as follows: _d_1 <-- _d_1' = _d_2/_u _d_2 <-- _d_2' = _d_1/_u where _u = det(_H) = 1 + (_d_1*_b_1*_b_1)/(_d_2*_b_2*_b_2) . The output value of _b is as follows: _1 _b_1 <-- _b_1' = _b_2*_u If rescaling is needed, SSRROOTTMMGG/DDRROOTTMMGG will further modify these output values before the end of the routine. See case 4. CCaassee 44:: Rescaling If the scaling factors become either very large or very small, the scaling and rotation operations may lose a lot of accuracy; therefore, each scaling factor from case 2 or case 3 is kept within the range: _g_a_m_m_a**(-2) <= | _d(_i)' | <= _g_a_m_m_a**2 , for _i = 1, 2 where _d(1)' = _d_1', _d(2)' = _d_2', and _g_a_m_m_a = 2.0**12 = 4096.0. At the end of case 2 or case 3 assignments, if either of the scaling factors falls outside this range, this routine must rescale that factor (and the corresponding elements of _H) to bring its size back within the specified range. 48 - log2(_g_a_m_m_a) = 48 - 12 = 36 bits (~ 10 decimal digits) Rescaling is performed as follows: If either _d_1 or _d_2 is 0, no rescaling is done; otherwise, let _q(_i) = int(log_base__g_a_m_m_a(sqrt(|_d(_i)'|))) = int(log2(|_d(_i)'|)/24) , for _i = 1 , 2. Then the following is true: _q(_i) < 0 , if |_d(_i)'| < _g_a_m_m_a**2 _q(_i) = 0 , if _g_a_m_m_a**(-2) <= |_d(_i)'| <= _g_a_m_m_a**2 _q(_i) > 0 , if |_d(_i)'| > _g_a_m_m_a**2 Furthermore, |_q(_i)| represents the number of times _d(_i)' must be multiplied (or divided) by _g_a_m_m_a**2 to return it to the proper range of values. In this case, the _r_p_a_r_a_m values set on output are as follows: _r_p_a_r_a_m(1) <-- -1.0 _r_p_a_r_a_m(2) <-- _h_1,_1' = _h_1,_1*_g_a_m_m_a**_q(1) _r_p_a_r_a_m(3) <-- _h_2,_1' = _h_2,_1*_g_a_m_m_a**_q(2) _r_p_a_r_a_m(4) <-- _h_1,_2' = _h_1,_2*_g_a_m_m_a**_q(1) _r_p_a_r_a_m(5) <-- _h_2,_2' = _h_2,_2*_g_a_m_m_a**_q(2) The output values of the scaling factors are as follows: _d_1 <-- _d_1'' = _d_1'*_g_a_m_m_a**(-2*_q(1)) _d_2 <-- _d_2'' = _d_2'*_g_a_m_m_a**(-2*_q(2)) The output value of _b_1 is as follows: _b_1 <-- _b_1'' = _b_1'*_g_a_m_m_a**_q(1) SSEEEE AALLSSOO RROOTT(3F), RROOTTGG(3F), RROOTTMM(3F) Gentleman, W. M., "Least Squares Computations by Givens Transformations Without Square Roots," _J_o_u_r_n_a_l _o_f _t_h_e _I_n_s_t_i_t_u_t_e _f_o_r _M_a_t_h_e_m_a_t_i_c_a_l _A_p_p_l_i_c_a_t_i_o_n_s 12 (1973), pp. 329 - 336. Lawson, C., Hanson, R., Kincaid, D., and Krogh, F., "Basic Linear Algebra Subprograms for Fortran Usage," _A_C_M _T_r_a_n_s_a_c_t_i_o_n_s _o_n _M_a_t_h_e_m_a_t_i_c_a_l _S_o_f_t_w_a_r_e, 5 (1979), pp. 308 - 325. This man page is available only online.